The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff chemical reaction network modeling the BCR signaling network from Barua et al.. We use ReactionNetworkImporters to load the BioNetGen model files as a Catalyst model, and then use ModelingToolkit to convert the Catalyst network model to ODEs.
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit gr() datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 # generate ModelingToolkit ODEs @timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) rn = prnbng.rn @timeit to "Create ODESys" osys = convert(ODESystem, rn) u₀ = prnbng.u₀ p = prnbng.p tspan = (0.,tf) @timeit to "ODEProb No Jac" oprob = ODEProblem(osys, u₀, tspan, p) @timeit to "ODEProb DenseJac" densejacprob = ODEProblem(osys, u₀, tspan, p, jac=true)
Parsing parameters...done
Adding parameters...done
Parsing species...done
Adding species...done
Creating ModelingToolkit versions of species and parameters...done
Parsing and adding reactions...done
Parsing groups...done
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 1122-element Vector{Float64}:
299717.8348854
47149.15480798
46979.01102231
290771.2428252
299980.7396749
300000.0
141.3151575495
0.1256496403614
0.4048783555301
140.8052338618
⋮
1.005585387399e-24
6.724953378237e-17
3.395560698281e-16
1.787990228838e-5
8.761844379939e-13
0.0002517949074779
0.0005539124513976
2.281251822741e-14
1.78232055967e-8
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, u₀, tspan, p, jac=true, sparse=true) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 359s / 99.2% 56.8GiB / 99.5%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
ODEProb DenseJac 1 298s 83.4% 298s 45.1GiB 79.9% 45.1G
iB
ODEProb No Jac 1 28.3s 7.93% 28.3s 5.54GiB 9.80% 5.54G
iB
ODEProb SparseJac 1 20.3s 5.70% 20.3s 4.21GiB 7.46% 4.21G
iB
Parse Network 1 7.25s 2.03% 7.25s 855MiB 1.48% 855M
iB
Create ODESys 1 3.22s 0.90% 3.22s 794MiB 1.37% 794M
iB
──────────────────────────────────────────────────────────────────────────
──
@show numspecies(rn) # Number of ODEs @show numreactions(rn) # Apprx. number of terms in the ODE @show numparams(rn) # Number of Parameters
numspecies(rn) = 1122 numreactions(rn) = 24388 numparams(rn) = 128 128
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
u = copy(u₀) du = similar(u) @timeit to "ODERHS Eval1" oprob.f(du,u,p,0.) @timeit to "ODERHS Eval2" oprob.f(du,u,p,0.) # force compilation for dense and sparse problem rhs densejacprob.f(du,u,p,0.) sparsejacprob.f(du,u,p,0.) J = zeros(length(u),length(u)) @timeit to "DenseJac Eval1" densejacprob.f.jac(J,u,p,0.) @timeit to "DenseJac Eval2" densejacprob.f.jac(J,u,p,0.)
ERROR: syntax: expression too large
Js = similar(sparsejacprob.f.jac_prototype) @timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.) show(to)
ERROR: syntax: expression too large
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol, legend=false, fmt=:png)